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Abstract. Including derivative information in the modelling of moving interfaces has been 
proposed as one method to increase the accuracy of numerical schemes with minimal additional cost. 
Here a new level set reinitialization technique using the fast marching method is presented. This 
augmented fast marching method will calculate the signed distance function and up to the second- 
order derivatives of the signed distance function for arbitrary interfaces. In addition to enforcing 
the condition ||V(/)||^ = 1, where (f) is the level set function, the method ensures that V(||V(/)||)^ = 
and VV(||V(/)||)^ = are also satisfied. Results indicate that for both two- and three-dimensional 
interfaces the resulting level set and curvature field are smooth even for coarse grids. Convergence 
results show that using first-order upwind derivatives and the augmented fast marching method result 
in a second-order accurate level set and gradient field and a first-order accurate curvature field. 
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1. Introduction. The fast marching method (FMM) was introduced by Sethian 
as an efficient method to solve general front propagation problems where the prop- 
agation speed is monotonic. Since its introduction the fast marching method has been 
successfully utilized in seismology [2 , photolithography [3] , medical imaging [11 [5] |6] , 
and as a component in other numerical schemes [71 18] . The fast marching method is 
also an extremely efficient way to compute the distance to an interface [9 . It is this 
last application which is the focus of this work. 

Recent attention has been focused on increasing the accuracy of the level set 
method by including level set gradient information [10 . Results for situations where 
the velocity field does not depend on the current interface show that the accuracy 
can be increased with minimal additional computational effort. Issues do arise when 
the velocity field depends on the current interface description. Take for example the 
modelling of vesicles in external fluid flows. The vesicle membrane exerts bending 
forces on the surrounding fluid [71 [Til ttl] • These forces are related to the curvature 
and the variation of the curvature along the vesicle membrane. Mathematically this 
means that the bending forces are fourth order derivatives of the level set function. 
The level set function should be as smooth as possible for all time and still be able to 
describe small scale features of the vesicle. This can be accomplished by periodically 
reinitializing the level set. As of yet a numerical method to accurately reinitialize the 
level set function and its gradients has not yet been developed. 

In this work the reinitialization of a arbitrary level set function is considered. In 
addition to obtaining the signed distance function of an interface the method will allow 
for the accurate calculation of up to second order derivatives of the signed distance 
function. This will result in the calculation of smooth curvature fields, which aids 
in the stability of numerical methods depending on this quantity and its derivatives. 
In Sec. [2] the original fast marching method for reinitialization is briefly presented. 
The augmented fast marching method is shown in Sec. [3] Here both the two- and 
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three-dimensional systems are considered. Two-dimensional convergence studies and 
results are presented in Sec. |4] while the three-dimensional results are shown in Sec. 

El 

2. The Fast Marching Method for Reinitialization. Consider an interface 
implicitly defined as the zero of a function, r{t) = {x : (l){x^t) = 0}, moving with a 
monotonic speed of F. The time, r, at which the interface crosses a point x is the 
solution to the Eikonal equation, F||Vr(ic)|| = 1. If the speed of the front is one {i.e. 
F = 1) then the time at which the interface will cross a point x is the distance from 
the point to the interface. By setting (j){x) = t{x) and solving ||V0(x)|| = 1 using the 
fast marching method it is possible to obtain the distance function of the interface. 
Denoting the region enclosed by the interface as the negative of the distance function 
results in a signed distance function. In this section the application of the original 
fast marching method to level set reinitialization is presented. More information can 
be found in references [l] |9] . 

2.1. Basics of the FMM. A two-dimensional computational domain with a 
uniform grid spacing of h has an embedded interface, P. The goal is to calculate 
the signed distance function to the interface without any spurious motion of P. To 
accomplish this the Eikonal equation with F = 1, 

\wn = 1, (2.1) 

is solved in the entire computational domain. In the fast marching method upwind 
derivatives are used to approximate V^. This enforces a causality on the propagation 
of information. Consider two points, x and y, where the point x is closer to the 
interface than point y. Due to the fact that F = 1 and the use of upwind derivatives 
the value at any given point only depends on those points closer to the interface. 
Thus the value (j){y) may depend on the value of ^(x), but the value of (/)(a?) will 
never depend on the value of (piy). This leads to an ordering of the nodes which 
needs to be maintained throughout the reinitialization procedure. 

When using the fast marching method three sets of nodes are maintained. The 
first are accepted nodes, A, which are those nodes where a value of (j) has already been 
calculated and accepted. The second set are trial nodes, T. The trial set contains 
those nodes which might next be added to the accepted set. Finally, the distant set, 
are those nodes which are too far from the interface to be added to the accepted 
set. Due to the causality of the FMM we are guaranteed that if a? G A, y G T, and 



z e D then |0(ic)| < < |^(^)|- See Fig. 2.1 for the relationship between the 

three set of nodes. 



• Accepted 
O Trial 
O Distant 




Fig. 2.1. The Accepted, Trial, and Distant set of grid nodes. 
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2.2. Initializing the FMM. To initialize the FMM mark all nodes in the do- 
main as in the distant set. Define a grid cell as the region enclosed by the four 
points Xij^ Xi-^ij, Xij-^i, and Xi-^ij-^i. Identify all grid cells which contain the 
interface and mark the nodes associated with these cells as in the accepted list. 

The standard fast marching method can be initialized by explicitly calculating 
the level set function for all nodes associated with cells containing the interface. De- 
fine a bicubic interpolation function, P(x), approximating the level set function in 
a given cell. In the absence of gradient information the bicubic function can be ob- 
tained by ensuring that P{Xm,n) = (t>m,n, dxPiXm^n) = dx(l)m,n, dyP{Xm,n) = dy(t)m,n, 

and dxyP{xm,n) = dxy4>m,n for m = z + 1 and n = j, j + 1 at grid cell ^i^j- All 
derivatives of the level set can be obtained by using standard finite difference ap- 
proximations. At every grid point Xi^j in the initially accepted the list the point y 
is calculated such that p{y) = and Vp{y) x {xi^j — y) = 0. The distance to the 
interface is then \\x — y\\2. This results in a second-order approximation to the true 
distance function [9]. 

2.3. Updating Nodes in the FMM. After the initial accepted list is deter- 
mined the remaining grid nodes are updated in an ordered manner. All grid nodes 
adjacent to the initial accepted list are given estimates of the distance function by 



solving an upwind discretization of Eq. (2.1), 



iDt4>)'+{D^<pf = l, (2.2) 

where and are one-sided derivatives. The appropriate derivative is chosen 
based on the direction of neighboring accepted nodes. Let the node being updated be 
Xij with nodes Xi-ij and Xij-^i in the accepted list. To first order the derivatives in 
this case would be D~(f) = — (j)i_i^j)/h and D^(j) = — (j)ij)/h^ where h is 



the grid spacing. In general solving Eq. (2.2) will result in two real roots, 0i and 



[9]. The smallest root which is larger than the surrounding stencil nodes is taken to 
be the accepted value. In the example case given the solution would be the smaller 
of 4>i and 4>2 that is larger than both (j^i-ij and 

The remaining grid nodes are updated in the following fashion. The grid node 
in the trial list with the smallest level set value is moved into the accepted list. All 
nodes surrounding the newly accepted node which are in either the trial or distant 



lists are updated by solving Eq. (2.2). Any updated nodes in the distant list are 
moved into the trial list. This procedure is repeated until no nodes remain in the trial 
list. The next node to be added to the accepted list is easily obtained if the trial list 
is maintained as a sorted list such as a heap. 

3. The Augmented Fast Marching Method. To improve the accuracy of 
fast marching based reinitialization schemes it is proposed to solve an extension of 



Eq. (2.1 ). Begin by taking up to second order derivatives of the square of the Eikonal 



equation with F = 1: 



V^- V0= 1, (3.1) 
V(V^- V^) = 0, (3.2) 
VV (V^ • V(^) = 0. (3.3) 
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In two dimensions this results in six equations, 









= 1, 


(3.4) 
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(3.6) 
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= 0, 


(3.9) 



with (px denoting partial derivative of (j) with respect to x. The values of interest are 
the level set function and up to the second derivatives of the level set function: ^, 
(j)x^ (l)y^ (j^xxi 4^yy^ and (j)xy. By using finite difference approximations to first order 
derivatives a set of six equations is obtained for the six unknowns. The particular 
discretization used in this work is given in Sec. |3.2[ while the extension to three 
dimensions is presented in Sec. |3.3[ 



3.1. Initializing the AFMM. Initialization of the augmented fast marching 
method proceeds in a manner similar to the standard fast marching method. Let the 
grid cell contain the interface. It is possible to define a bicubic interpolant over 
this cell using the given data. As this work was designed to work with the gradient 
augmented level set method it is assumed that the level set value, ^, and the gradient 
of the level set, = i/; = (?/;^,?/;^), is available at the four grid nodes associated 
with Qij. To compute the bicubic interpolant it is necessary to define (j)xy at the four 
grid points. This value is calculated as the average of the derivatives of the gradient 
field. At any point (j)xy = {Dxi^'^ -\- Dyip^)/2 where Dx and Dy represent the centered 
second order finite difference approximations to the first derivative. 

To initialize the augmented fast marching method all nodes associated with grid 
cells containing the interface are moved into the accepted list. Values for the level set, 
^, the gradient of the level set, -i/?, and the Hessian of the level set, are calculated 
at each of these initially accepted points. 

The initialization follows a technique developed for accurately calculating the 
curvature of level set functions [13 . Let the grid node Xij be in the initially accepted 
list. A 3 X 3 sub-grid is centered at Xij. The spacing of this sub grid is taken to 



be ah^ where a < 1 and h the uniform grid spacing, see Fig. 3.1 Each of the nine 
points in the sub-grid have a signed distance function value calculated by minimizing 
the function 111/ — iCo||2 subject to P{y) = 0, where P{x) is the bicubic interpolant 
over the grid cell Qij and Xq is the sub-grid point. The required derivatives are then 
obtained by standard second-order finite difference schemes using the sub-grid data. 
Due to the small sub-grid size the accuracy of the initialized values is extremely high, 
see the results sections for details about the convergence rate. 

3.2. Updating Nodes in the AFMM. The remaining grid nodes are updated 
in an ordered manner. Nodes adjacent to the initially accepted list are updated by 



solving the discretization of Eqs. (3.1 )-(3.3). These nodes are placed into the trial list 



which is kept as a heap sort to ensure the fast retrial of the node with the smallest 
level set value. The node in the trial list with the smallest value is placed into the 
accepted list and any non- accepted adjacent nodes are updated. If an adjacent node 
is in the distant list it is moved into the trial list. This procedure is repeated until no 
nodes remain in the trial list. 



Fig. 3.1. The initialization grid at a grid point Xij with a sub-grid centered at the node. The 
distance from each node of the sub-grid to the interface is calculated and shown for three nodes. The 
values of (j), ip and H are then computed using finite difference approximations on the sub-grid. 



The updating procedure for the AFMM is more involved than the standard FMM. 
In the standard FMM a single nonlinear equation needs to be solved, Eq. (2.2). In 



the 2D AFMM there are six nonlinear equations. The solution obtained will depend 
on the particular discretization of the equations given in Sec. [3) In the case of two 
accepted nodes, one in each Cartesian direction, the discretization chosen is 
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(3.10) 










= 0, 


(3.11) 










= 0, 


(3.12) 




h H^yR^y - 






= 0, 


(3.13) 


ffyyjjyy _ 


f H^yR^y - 






= 0, 


(3.14) 


j^xx j^xy 


\- H^yu^y - 




i>y{D^H''y) 


= 0. 


(3.15) 



The operator and represent the appropriate one-sided derivatives at the grid 



point to be updated, see Sec. 2.3 or Ref. [9 



Also note that due to symmetry 
there are three components to the Hessian matrix which we denote as (j)xx = H^^^ 
(j)yy = Ryy ^ and (j)xy = H^y . It is worth noting that this particular discretization 
allows for the solution of the level set and gradient field first, Eqs. 



(3.10)-(3.12). 



Once a valid solution set {^,1/?} is calculated the solution of the Hessian field, Eqs. 
(3.13)-(3.15) can be determined. 



To determine if a calculated solution to the gradient system, Eqs. (3.10)-(3.12), 
is valid acceptance criteria similar to the classical FMM must be implemented. Due 
to the larger amount of information the acceptance criteria of the classical FMM is 
augmented. In addition to requiring that the solution level set value be larger than 
either of the grid point's neighbors it is also required that the calculated gradient be 
in the same general direction as the neighboring nodes. For example, let the nodes 
Xi-^ij and Xij-^i be accepted nodes when updating the value at Xij. In this case the 
chosen finite difference approximations would be the stencils for the positive one-sided 
derivatives in each direction, and . Let and tj^ = (^tjj^^tjjy^ be solutions 
to Eqs. (3.10)-(3.12 ). A valid solution set would satisfy <p > (j) > 
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• ^ and -i/? • > 0. If the solution set 1/?} satisfies these conditions 

then it is taken as vahd. 

When updating the values at Xij it is possible that only a single neighboring 
node is in the accepted list. In this case the discretization is modified to account for 
the reduced amount of information. To use only information in the x-direction the 
system becomes 





= 1, 
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(3.21) 



The general idea is to replace derivatives in the direction not present (^/-direction in 
this case) with the equivalent x-direction derivatives. A similar formulation can be 
made for using information only from the ^-direction and is given in Appendix [A) 

3.3. Extension of the AFMM to 3D. The extension to three dimensions is 



straightforward. Explicitly writing Eqs. (3.1)-(3.3) results in ten equations 
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(3.31) 



The quantities of interest are the level set function and derivatives of (j) up to 
second order. Using upwind finite difference approximations for the first derivative 
this results in a set of ten nonlinear equations for the ten unknowns. 

In the three dimensional case initialization of a grid point Xij^k occurs over a 
3x3x3 sub-grid centered on Xij^k- To use numerical approximations of the derivatives 
it is not necessary to determine the level set value for all 27 points on this sub- 
grid as only 19 points are used during the derivative calculations. As in the two- 
dimensional case these sub-grid points have their distance to the interface calculated 
by minimizing \\y — XqW^ subject to P{y) = 0, where P{x) is the tricubic interpolant 
over the grid cell ^ij^k- When computing the tricubic coefficients it is sufficient to 
ensure that the following values are satisfied at the eight corner nodes of the grid 
cell: 0, (j)x^ (j)y^ (j)zi (j^xyi 4^xzi 4^yzi ciud (pxyz- Assumiug that only the level set and 
the gradient of the level set are available during reinitialization any higher order 
derivatives are obtained by averaging the appropriate derivatives. For example the 
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third-order derivative necessary for the tricubic function is given by (j)xyz = {Dyz^p^ + 
Dxzi^'^ + Dxyi^^)/^^ where D^y is the finite difference approximation to dxy 

Updating the remaining grid points proceeds in the same fashion as the two- 
dimensional case. Due to the higher number of dimensions there are seven possibil- 
ities from where information will be propagating. The specific discretizations for all 
possibilities have been presented in Appendix [B) 

3.4. The AFMM Algorithm. One consequence of choosing the above dis- 
cretizations is that the solution for the gradient system, Eqs. (3.1Q)-( [3?T2 ) for the two 
dimensional case, can be computed first. Once a valid solution set {^, i/^} is calculated 
the solution for the Hessian field, Eqs. (3.13)- (3. 15) can be determined. In practice 



this is done by first computing for the gradient field in the entire domain. This not 
only results in level set and gradient functions at every grid point, it also determines 
the proper ordering to ensure that upwind derivatives are calculated correctly. The 



Hessian field is then obtained by solving Eqs. (3.13)-(3.15) using this pre-determined 
ordering. 

The algorithm for the augmented fast marching method can be summarized 
thusly: 

1. Mark all nodes as in the Distant list. 

2. Initialize all nodes associated with cells containing the interface by explic- 
itly solving for the signed distance function on a sub-grid centered at the 
node. Derivatives of the level set are calculated by standard finite difference 
approximations on this subgrid. Move these nodes into the Accepted list. 

3. For all nodes in the Distant list that lie next to a node in the Accepted list 



calculated updated values by solving Eqs. (3.4)-(3.6) in 2D or Eqs. (3.22) 



(3.25) in 3D. Move these nodes into the Trial list. 



Select the node with the smallest level set value in the Trial list and move 
it into the Accepted list. All nodes next to the newly accepted grid point in 
either the Trial or Distant list have updated values calculated by solving Eqs. 



(|3^-(|3^ in 2D or Eqs. ( |3.22D -( |3.25D in 3D. Any node that is updated and 
in the Distant list is moved into the Trial list. 
Repeat Step 4 until the Trial list is empty. 

Calculate the Hessian field by solving Eqs. (3.7)-(3.9) in 2D or Eqs. (3.26)- 



(3.31) in 3D. Update nodes using the same order as they were added to the 



Accepted list. 



4. Two Dimensional Results. Here convergence and sample results are pre- 
sented for two dimensional interfaces. All numerical derivatives were calculated using 
standard first-order one-sided finite difference schemes. The domain is the region 
[—2,2]^ and a uniform grid spacing h is used. No domain boundary conditions are 
needed as all information flows from the interface outwards. In all cases the the inter- 
face is initially described by a level set function 0o and its gradient fleld, = V^o- 
The sub-grid is taken to have a spacing of O.lh. The SLSQP algorithm |T4l [151 of 
the NLopt software library [16 was used to determine the closest point on the inter- 
face during initialization. All nonlinear systems were solved using the GNU Scientific 
Library [17 . 

4.1. Investigation of the Expected Errors for the Gradient System. 

The accuracy of a fast marching method depends on the order in which nodes are 
updated. A standard error analysis is difficult due to the nonlinear nature of the 
systems involved. Instead sample analytic solutions for a circular interface of radius 
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Ro are considered. This interface is to be reinitialized on a grid with a uniform spacing 
of h. Consider the AFMM updating of a grid point given by Xij = (x, x), x > 0, and 
V2x^ < Rq. This point hes on the domain diagonal. In this situation it is known 
that the grid points to the right, and above, Xij-^i^ are closer to the interface 

than Xi^j, Fig. |4.1 If the neighboring nodes have the exact solution prescribed the 
can be investigated. 



Fig. 

error of node x 




x^: Node to be update 

Xi+i J and Xi j+i: Nodes with exact solution 



Fig. 4.1. Sample stencil used for analysis of expected errors. The grid point 



updated by the AFMM using the nodes Xi-^ij and Xij-^i. 
exact solution prescribed. 



The nodes Xi-^ij and Xij-^i have the 



Consider the gradient system, Eqs. (3.4)-(3.6). Using the exact solution for nodes 



Xi-^ij and Xij-^i the level set, (/), and gradient vector, V0 = = [ip^^ipy)^ at Xi 



can be calculated as 



The true solutions are (ptme 
of the level set is calculated as (j) 



2xVh^ + 2hx + 2x2 
" h^2x 
h^2x 



^0, 



(4.1) 
(4.2) 



Ro and V^f^^^ = V^J^^^ = The error 



Rq while the error of the gradie nt v ector 
where (/), ip^ 



(4.2). These errors are shown in Fig. 4.2 for locations ranging from x 



and ipy are given in Eqs. (4.1) and 
to X = 0.5 



and grid spacings of h = 0.1, 0.05, and 0.025. 
From the results shown in Fig. 



4.2 it becomes apparent that as the grid size 



decreases the overall error decreases for both the level set and gradient vector. The 
maximum error of the level set function observes a first-order convergence. The gra- 
dient vector, on the other hand, has a fixed error of {V2-1) /V2 as X approaches 
zero irregardless of the grid spacing. This result should not be unexpected. As x 
approaches the origin the variation of the gradient field increases inversely to the dis- 
tance from the origin. For example, given a circular interface and a point x = (x, x) 
the exact variation of the x component of the gradient field in the x-direction is given 
dxi^^ = 1/ (2V^x). As the distance from the origin decreases due to a smaller grid 
spacing the variation in the gradient field will increase by an inverse amount. This 
results in a fixed error being introduced into the system near the origin. Despite this 
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0.30f 




(a) Level Set Error (b) Gradient Error 



Fig. 4.2. The error of the level set and gradient functions based on the stencil given in Fig. 
\4.l\ f0r grid spacings of h = 0.1, 0.05, and 0.025. 



0(1) error of the gradient field the overall error decreases rapidly. As the Hessian 
system depends on the solution of the gradient vector it should be expected that the 
errors for the Hessian field, and thus the curvature, will increase as one approaches 
the origin. 

To conclude this section it should be noted that the error at any given point will 
decrease for both the level set and gradient vector. As an example the error for the 
level set and gradient vector are shown in Fig. |4.3| for the points x = 0.01, 0.1, and 0.2 
using grid spacing ranging from 10~^ to 10~^. For the three points considered second- 
order convergence is observed for both the level set and gradient vector solutions. 




Grid Spacing 
(a) Level Set Error 



Grid Spacing 

(b) Gradient Error 



Fig. 4.3. The error of the level set and gradient functions based on the stencil given in Fig. 
\4 A cit the points x = 0.01, 0.1, and 0.2 for grid spacings ranging from 10~^ to 10~^. 



4.2. Accuracy of the Initialization Method. The accuracy of the initializa- 
tion scheme presented in Sec 13.11 is shown here. A circle of radius 1 with an initial 
level set given by = — e is considered. The resulting convergence rate is 

seen in Fig. |4.4| 

It is observed that the convergence rate for the level set is approximately 4^^- 
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order, for the gradient field it is 3^^ order, for (jy^x and (j)yy 2^^-order convergence is 
seen, and the (j)xy value is 1^^-order accurate. 

4.3. Accuracy of the 2D AFMM. First consider a unit circle with the initial 
level set of (/)o = +^ — e. Sample level set and curvature results at grid sizes ranging 



from 20^ to 500^ are shown in Figs. 4.5 and 4.6 The level set contours shown utilize 
the additional information provided by knowledge of the derivatives of the level set 
function. Even in the coarsest mesh, 20^, the level set is extremely smooth and the 
curvature field is smooth outside of the interface. At such a coarse mesh there are 
not enough grid points to accuractly describe the large variations of curvature which 
occurs in the region given by (/> < 0. As the number of grid points increases this error 
vanishes. 



Considering the results shown in Sec. 4.1 errors will be reported for both the 
entire domain and for the local region around the interface given by |^| < 9h. This 
width was chosen based on the stencils needed for a 5^^-order WENO local level set 
scheme [18 . To account for the fact that the curvature grows as Xj \J ^ as one 
approaches the center of the given level set function all curvature errors are reported 
d.^ {n — Ke) / whcrc K is the curvature using the AFMM and is the exact curvature. 



The errors for the unit circle are given in Figs. 4.7 and 4.8 

Despite the use of first-order finite difference approximations the resulting level 
set is second-order accurate and the gradient is 0{h^^'^) in the entire domain while 
the curvature is 0{h). When only considering the region close to the interface this 
increases to 0{h^^'^) for both the level set and gradient fields and 0{h^^'^) for the 
curvature field. 



The convergence results verify the conclusion obtained in Sec. 4.1 The large 
variation in the gradient vector as one approaches the center of the circle creates a 
uniform error that can not be overcome. The fact that the L2-norm does converge 
indicates that this error is localized around the origin. This is further supported by 
the convergence of both the L2 and L^o-norms in the region given by < 9h. 
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-2-1 1 2-2-1 1 2 

(a) 202 (b) 502 




-2-1 1 2-2-1 1 2 

(c) 2002 (d) 5002 



Fig. 4.5. The level set after reinitialization for a unit circle with an initial level set of (f)o = 

2_|_ 2 

-\-y _ g_ ^'/le interface is given by the thick red contour. These contours utilize the sub-grid 
information provided by the additional derivative information. 

Next consider an elliptical interface with a major axis of 1.5 and a minor axis of 
0.5: (j)o = (x/1.5)^ + {y/0.5)'^ — 1. Results for various grid sizes are presented in Figs. 
14.91 and [4. 101 As for the unit circle the level set field is smooth even at coarse meshes. 
As the grid becomes finer the level set field becomes smoother, particularly along the 
long axis of the ellipse. The curvature field is also smooth, even for the extremely 
coarse grid. For all of the grids the curvature field has the same qualitative shape. 
An issue with the curvature can be observed towards the center of the ellipse. This 
is due to the gradient vector field having a sharp discontinuity in that region. This 
is illustrated in Fig. |4.11[ which shows the right half of the reinitialized solution on 
a 50^ grid. The resulting gradient vector field switches sign when crossing the x-axis, 
resulting in errors being introduced into the Hessian, and therefore the curvature. 
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solution. 

As in the unit circle case the convergence of the elliptical interface is presented 
for both the entire domain and the region given by |0| < 9h. The exact solution was 
calculated by explicitly determining the signed distance function on a 4000^ grid. The 
convergence results for the elliptical case match those of the unit circle. 

This section concludes by presenting three additional sample interface. In each 
case a coarse grid of 50^ is compared to a fine grid of 500^. Results are for two circles 
of radius 0.75 centered at (0.8125,0.4125) and (-0.8125,-0.4125), a Cassini oval 
with an initial level set of = {{x - a)^ + y'^) + ((x + a)^ + y'^) - with a = 0.99 
and b = 1.01, and a star interface given by = x'^ -\- y'^ — 1 + sin(5^)/4 with 



ArcTan(?//x), see Figs. 4.14 to 4.16 
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10'^ 10'^ 10'^ 10" 

Grid Spacing Grid Spacing 

(a) L2 Error (b) Loo Error 



Fig. 4.7. The error in the entire domain for the unit circle. The Loo error does not converge 
for the gradient and curvature fields, as expected by the results of Sec. \4 A 




10'^ 10'^ 10'^ 10" 

Grid Spacing Grid Spacing 

(a) L2 Error (b) Loo Error 



Fig. 4.8. The error for the unit circle in the region given by < 9h, where h is the grid 
spacing. In this case both the L2 and Loo error-norms converge. 



In all three cases level set and curvature fields at the coarse grid qualitatively 
match those at the finer grid. Errors in the curvature field are again demonstrated 
in regions where level set fronts collide, such as the diagonal in the two circle case of 
Fig. 4.14[ the y-axis of the Cassini oval, Fig. 4.15[ and along several of the diagonals 
of the star shape. Fig. 4.16 These error diminish as the grid is refined. 



5. Three Dimensional Results. Sample three-dimensional results are pre- 
sented here. Due to computational limitations only a limited set of convergence 
results will be presented for three-dimensional surfaces. The domain will be the cube 
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(d) 5002 



Fig. 4.9. The level set after reinitialization for an elliptical interface with an initial level set of 
(f)o = (a;/1.5)2 + (y/0.5)2 — 1. The interface is given by the thick red contour. These contours utilize 
the sub-grid information provided by the additional derivative information. 



of [—2, 2]^ with a uniform grid spacing of h in all directions. 

First consider a spherical surface with a radius of one. The resulting isosurfaces 
for grids of 25^ and 100^ are presented in Figs. 



5.1 



and 5.2, respectively. The AFMM 



on both grids results in a smooth level set and mean curvature field. 



and 5.4 



Limited convergence results for the spherical interface are presented in Figs. 5.3 
As seen in the two-dimensional results the L 



for the gradient and curvature fields, as expected from Sec. 4.1 



norm errors do not converge 
It appears that 



in the L2 error for the entire domain converges at third-order. This is most likely 
due to the additional directions that information may travel. The most error will be 
introduced when only a single neighbor node is available during the calculation of 
updated values. This is much less likely to occur in three dimensions than in two, 
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Fig. 4.10. 

level set of (f)o 
contour. 



The curvature field afl:er reinitialization for an elliptical interface with an initial 
- (x/1.5)^ + (y/O.S)^ — 1. The curvature of the interface is given by the thick red 



resulting in an increase in the accuracy of the scheme. 

The result for an ellipsoidal surface give n by = \/(x/1.6)^ + (^/1.2)^ + (z/0.5)^- 
1.0 using a 100^ grid is shown in Fig. 5.5 while that for a three-dimensional Cassini 
oval given by 0o = {{x - a)^ + + z^) ({x + a)^ + + - with a = 1.29 and 
b = 1.3 on a 100^ grid is shown in Fig. 5.6 In both cases the level set and overall 



curvature field are smooth. As in the two-dimensional cases slight issues are observed 
for the Cassini oval in regions where level set contours collide. 

6. Conclusion. In this article the augmented fast marching method for reini- 
tialization of level sets is presented. This work builds upon the fast marching method 
work of Chopp [9] and the gradient augmented level set work of Nave et. al [10 . This 
method increases the accuracy of standard level set schemes by calculating the signed 
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Fig. 4.11. Vector field oj the reinitialized elliptical interface solution on a 50^ grid. A sharp 
discontinuity of the vector field occurs along the x-axis. 
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10 

Grid Spacing 
(a) L2 Error 



10 

Grid Spacing 
(b) Loo Error 



Fig. 4.12. The error in the entire domain for the ellipse. The Loo error does not converge for 
the gradient and curvature fields, as expected by the results of Sec. \4 A 



distance function and up to second-order derivatives of a general interface. Results 
show that both the level set and curvature fields are smooth for a wide variety of 
interfaces and in both two- and three-dimensions. It has also been demonstrated that 
the scheme calculates a higher than second-order accurate level set and gradient vec- 
tor field while the resulting curvature is slightly higher than first-order accurate. This 
was accomplished using standard first-order upwind derivatives, as in the original fast 
marching method. Unlike partial differential equation based reinitialization schemes 
a fast marching based method has the advantage of only requiring a single pass to 
update a level set far from the signed distance function. 

The additional accuracy has already proven useful in practice. This technique has 
been used in the modelling of vesicles in general fiows [71 [19] . In such a simulation 
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-2 -2 

10 10 
Grid Spacing Grid Spacing 

(a) L2 Error (b) Loo Error 



Fig. 4.13. The error for the ellipse in the region given by |^| < 9h, where h is the grid spacing. 
In this case both the L2 and Loo error-norms converge. 

second-order derivatives of the curvature are required. The use of this reinitiahzation 
technique, in conjunction with the gradient- augmented level set method, allowed for 
meaningful simulations to be performed using a single workstation. 
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Fig. 4.14. Contours of the level set and curvature for dual circles centered at (0.8125,0.4125) 
and (—0.8125, —0.4125) and both with radius of 0.75 on a coarse and fine grid. 




Fig. 4.15. Contours of the level set and curvature for a cassini oval with an initial level set of 
(f>0 = {(x — a)^ + y-^) ((x + a)^ + y^) — 6^ with a = 0.99 and b = 1.01 on a coarse and fine grid. 




Fig. 4.16. A star interface given by (j)o = \J ^ — 1 + sin(5^)/4 with = ArcTan(y/x) on 
a coarse and fine grid. 
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(a) Level Set (b) Curvature 

Fig. 5.2. Sphere of radius one on a 100^ grid. The red isosurface represents the value at the 
interface. 
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0.02 0.04 0.06 0.1 0.02 0.04 0.06 0.1 

Grid Spacing Grid Spacing 

(a) L2 Error (b) Loo Error 



Fig. 5.4. The error for the spherical surface in the region given by \(f)\ < 9h, where h is the 
grid spacing. In this case both the L2 and Loo error-norms converge. 



(a) Level Set 



Fig. 5.5. Ellipsoid given by (f)o = 
red isosurface represents the interface. 



(b) Curvature 

^/{x/1.6y + (2//1.2)2 + (z/0.5)2 - 1.0 on a 100^ grid. The 




(a) Level Set (b) Curvature 

Fig. 5.6. Cassini oval given by (j)o = ((x — a)^ + 2/^ + z^) {(x + a)^ + 2/^ + z^) — "w^^^^ a = 
1.29 ano? 6 = 1.3 on a 100^ grid. The red isosurface represents the interface. 
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Appendix A. Two dimensional discretization using only y-direction in- 



formation. The two-dimensional discretization of Eqs. (3.4)- (3.9) using only infor- 
mation from the y-direction is 



1 



j^yy^yy ^ j^xyj^xy ^ ^^D^H'^y) + iljy{D^Hyy) = 

jjxxjjxy ^ jjyyjjxy ^ H'"'') + ^ljy{D^H''y) = 



Appendix B. Three dimensional discretizations. The full discretizations 
for all cases in three dimensions are included here. There are a total of seven cases to 
consider. 

B.l. Information from the x-, y-, and z-directions. 

riDt^) + riD^^) + riDf^) = 1 
riDtr) + r{D^r) + notr) = o 

riD^r) + r{D^r) + notr) = o 

jjxyjjxy jjyyjjyy _^ Jjyz jjyz _^ ^^(D±7jyy) + tf^V {D^ H^y) + Xp^D^ HVy) = 

^xz^xz ^ jjyzjjyz _^ jj^^jj^^ _^ iP''{D^H'') + x(}y{D^H'') + ip^DfH'') = 

jjxxjjxy _^ jjyyjjxy _^ jjxz jjyz _^ H'^V) + Ip^iD^H^'y) + tp^{DfH''y) = 

jjxxjjxz jjzzjjxz ^ jjxyjjyz _^ ^^(^D^ff^^) + tl^^ {D^ H'''') + V'^(Z)±iJ^^) = 

HWRy^ + H'^Ry' + H'''H''y + ip^'iDtHy') + tpy{D^Hy') + ip''{DfHy') = o 

B.2. Information from the x- and y-directions. 

riDtr) + r{D^r) + r{Dtr) = o 
riDtr) + r{D^r) + r^o^r) = o 

jjxxjjxx ^ jjxyjjxy _^ Jjxzjjxz _^ ^^(£,±iyxa;) _^ + iP^DtH^') = 

jjxyjjxy jjyyjjyy _^ + tl^'' {D^ HyV) + 1py{D^Hyy ) + %p\D^ Hy") = 

iJ^^iJ^^ + ijs'^ijf^ + iJ^^iJ^^ + ip''{DtH'') + il^y{D^H'') = 

ipy{D^H^y) + + {D^Hy'))/2 = o 
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B.3. Information from the x- and z-directions. 

jjxxjjxx ^ jjxyjjxy ^ jjxz jjxz ^ ^^(£)±jy^^) + ij,y{D^H''y) + (D^il^^) = 
jjxyjjxy jjyyjjyy _^ jjyzjjyz _^ ^^(£)±^TO) + = Q 

iP^'iD^H'^') + i^motH'^y) + (D±J?f"))/2 + iP'iDfH''') = 
^ra^^y^ + iJ^^iJf^ + iJ^^^iJ^y + xp='{D^Hy') + tpy{DfHyy) + ip^D^Hy^) = 

B.4. Information from the y- and z-directions. 

riDtr) + r{D^r) + notr) = o 

jjxx^xx ^ jjxyjjxy _^ ^xz^x^ _^ ^y (D^ H'^^) + tl^^ {Of H'"' ) = 
^xy^xy _^ ^ra^yy _^ JJVz jjyz _^ ,/,»^(£)±iJ^y ) + tf^V {D^ Ryy) + Xp'{DfHyy) = 

^xz^xz ^ ^j/z^!/z + il^^il^^ + t(,''{DfH''') + ipy{D^H'') + xp^DfH'') = 

jjxxjjxy ^ ^w^xy ^ j;^xz^yz ^ ^^(£,±^^^^) + (D^ H^'y) + lP'{DfH''y) = 
^xx^xz _^ ^zz^xz _^ ^xy^yz _^ ^^(£)±ij-^^) + xj^V {D^ H''') + = 

jjyyjjy^ ^ jjzzjjyz _|_ jj^^jj^v^ 

i^'^iiDfH''') + {DfH''y))/2 + ipy{D^Hy') + ti^'{DfHy') = o 

B.5. Information from the x-direction. 

(D^4>)7p'' + 'tpy'tpy + 'tp''tp' = 1 

t/;^(Z)±V^) = 
= 

jjxxjjxx ^ ^xy^xy ^ jjxzjjxz ^ ^^(£)±ij-^^) + ■^V {D^ H^'y) + tj}' {D^ H'^'') = 

^xy^xy _^ ^yy^yy _^ ^yz^yz _^ (£)±iJ-yy ) = Q 

^xz^xz _^ Hy^ny' + ff^^ff^^ + tp''{DtH'') = o 

jjxx^xy _^ jjyy^xy _^ jjxzjjyz _^ ^^(/jijj^y) + i^V {D^ HVy) + Ip^D^Hy') = 
^xx^xz _^ ^zz^xz _^ ^xy^yz _^ ^^(£,±ij^^) + tpy{D^Hy^) + l/^^(L>±iJ^^) = 
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B.6. Information from the y-direction. 
B.7. Information from the z-direction. 

^^^^ _^ + {Df(t))tl^' = 1 

jjxxjjxx ^ jjxyjjxy ^ jjxz jjxz ^ (^JJ^ H'"'') = 
ffxyjjxy jjyyjjyy _^ jjyzjjyz _^ ^^^^fffVy^^ = Q 

^xx^x^ + ^^^^x^ ^ ^^2/^2/^ ^ + %l^y{DfH''y) + tl^'iDfH''') = 

ij^^ij^^ + ij^^ij?/^ + _^ ^^^DfH'^y) + ^i^y^DfRyy) + ^i^'iDfRy) = o 
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